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The Hamiltonian Mean Field (HMF) model of coupled inertial, Hamiltonian rotors is a prototype 
for conservative dynamics in systems with long-range interactions. We consider the case where 
the interactions between the rotors are governed by a network described by a weighted adjacency 
matrix. By studying the linear stability of the incoherent state, we find that the transition to 
synchrony occurs at a coupling constant K inversely proportional to the largest eigenvalue of the 
adjacency matrix. We derive a closed system of equations for a set of local order parameters and 
use these equations to study the effect of network heterogeneity on the synchronization of the 
rotors. We find that for values of K just beyond the transition to synchronization the degree of 
synchronization is highly dependent on the network’s heterogeneity, but that for large values of 
K the degree of synchronization is robust to changes in the heterogeneity of the network’s degree 
distribution. Our results are illustrated with numerical simulations on Erdos-Renyi networks and 
networks with power-law degree distributions. 


I. INTRODUCTION 


The Hamiltonian Mean Field (HMF) model [1-6] is 
a paradigmatic model for conservative systems exhibit¬ 
ing long-range interactions. Examples of such systems 
include free electron lasers [7], rarehed plasmas [3], grav¬ 
itational n-body problems [8], etc. This model has at¬ 
tracted attention due to its striking dynamical properties 
which include second order phase transitions and violent 
relaxation towards persistent meta-equilibrium states [9] . 

The generalized HMF model describes the dynamics 
of N interacting rotors with phase angles and angular 
momenta {(d„,p„) : n = 1,2,... N} through the Hamil¬ 
tonian 
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Here the first sum represents the kinetic energy of the 
rotors with moments of inertia and the second the 
potential energy of coupling through a network adjacency 
matrix A where Anm 7^ 0 if there is an edge from node 
m to node n and Anm = 0 otherwise. For Anm > 0, the 
potential energy due to the interaction of rotors n and m 
is minimized when they are aligned, = 9m- Without 
loss of generality A can be taken to be symmetric, = 
A, since the asymmetric part of A does not contribute to 
the interaction term in (1). The overall coupling strength 
is represented by K and it is scaled by 1/N so that the 
energy per rotor has a finite limit as N ^ 00 . 

While the HMF model has been proposed as a model 
for systems with long-range interactions, in its commonly 
studied form these interactions are assumed to be of such 
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long range that the rotors are all-to-all coupled [Anm = 1 
in (1)]. A natural question is what is the effect of al¬ 
lowing a more general form for the interaction network. 
Such a generalization would include spatially distributed 
systems with decaying interactions, varying interaction 
strengths, and arbitrary interaction structure. Since the 
case of heterogenous moments of inertia was considered 
by two of the authors in [10], we will assume in the 
current paper that the moments of inertia are identical, 
7„ = 1, focusing on the effects of the network structure on 
the dynamics. Most previous studies for the HMF model 
except [8, 10-12] have considered the all-to-all case with 
Anm = 1 in (!)■ Chavanis et al. [8] consider stellar 
(gravitational) systems with interactions depending on 
the mass M = / of the stars and thus Anm = MnMm- 
Restrepo and Meiss [10] study the disordered HMF model 
where Anm = o-nO-m, and a and I have independent, het¬ 
erogeneous distributions. In terms of the network struc¬ 
ture, both of these variants of the HMF model can be 
thought of as dynamics on a weighted, all-to-all network. 
Ciani, Fanelli, and Ruffo [11] studied the HMF model on 
Erdds-Renyi networks. Another generative model, the 
Watt-Strogatz small-world network [13], was used by Ni- 
gris and Leoncini [12]. Both [12] and [14] obtain a de¬ 
scription of the dynamics in terms of the network model 
parameters that requires a model-fitting step. In contrast 
to these previous approaches that study specihc network 
ensembles, in this paper we will develop a more general 
theory that applies to any given network described by 
its adjacency matrix A. To test our theory, we will use 
Erdos-Renyi networks and networks with heterogenous 
degree distributions, such as networks with power-law de¬ 
gree distributions, i.e., “scale-free” networks [14]. While 
we will use both Erdos-Renyi and scale-free networks in 
our examples, we emphasize that our analysis does not 
rely on an assnmed generative mechanism for the net¬ 
work: it works directly with the network adjacency ma¬ 
trix. 

Below we determine onset of instability of the incoher- 
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ent state, obtain a self-consistent equation for a set of 
local order parameters and quantify the degree of syn¬ 
chrony in terms of a macroscopic global order parameter 
R. As in previous studies on network synchronization, 
e.g. [15], we find that the principal eigenvalue A of the 
network adjacency matrix is a key quantity in determin¬ 
ing the onset of synchronization. Finally, we quantify 
the maximum achievable synchrony for a given network 
structure and find that this maximum value is very robust 
to the heterogeneity of the network’s degree distribution. 

The rest of the paper is organized as follows. The 
model and its governing dynamical equations are de¬ 
scribed in §11. In §III we discuss the linear stability of 
the incoherent solution. We use this analysis to find the 
critical value of the coupling constant for the onset of 
synchronization. We then study the synchronized state 
in §IV and give results for the global order parameter as 
a function of the coupling strength in terms of a set of 
self-consistent equations for the local order parameters. 
We provide approximations to the solution of these equa¬ 
tions just past the onset of synchrony and in the strong 
coupling limit. Finally, we discuss our results in §V. 

II. NETWORK HMF MODEL 

In the original HMF model, and in most subsequent 
studies [1-4], all rotors in (1) were assumed to have the 
same moments of inertia, /„ = 1 , and the coupling was 
assumed to be all-to-all with equal strength, Amn = 1 - 
While such a simplified setting provides many insights, 
interactions are rarely uniform and all-to-all in practice. 
For example, the HMF model is a simplified model for 
an n-body gravitational system in one spatial dimension 
with periodic boundary conditions, keeping only one har¬ 
monic of the potential [ 6 , 8 ]; in this case, the interaction 
strength should be proportional to the product of the 
particle masses and decay with the separation of the par¬ 
ticles. 

With this motivation we allow for a general adjacency 
matrix, A, in (1), but simplify by setting /„ = 1. The 
resulting dynamical system is 

0n=Pn , (2) 

Pn — ^ ^ Anm sin {9m ^n) • ( 8 ) 

m—1 

As is usual, it is convenient to define order parameters 
to quantify synchronization. When the network is het¬ 
erogeneous, one can define a set of real, local order and 
phase parameters, {{Rn,'4’n) : n = 1 ,... N}, by 
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that characterize the coherence of inputs to a given node. 
Using these, (3) becomes 

Pn — KRn sici (V^n ^n) (5) 


The overall synchrony of rotors can be measured by a 
global order parameter [15] 
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Here || ... || denotes the average over nodes, 

n—l 

and dn denotes the effective degree of the node, 

N 

dn = ^ ^ -^nm ■ 
m—l 

The normalization in (6) is chosen so that i? = 1 if all 
rotors are in synchrony (0„ = 9m)- 

III. LINEAR STABILITY ANALYSIS 

In this section we study the incoherent state, in which 
the local order parameters are approximately zero 
and the rotors evolve approximately independently of 
each other, i.e., (2) and (3) become = Pn, Pn = 0. 
In this case = p„(0)< -I- 0„(O). Assuming that the 
initial momenta differ, p„(0) yf _Pm(0) for m ^ n, then 
each oscillator has a different frequency and (4) gives 
(l^np)i = Here (...)t denotes a time 

average, 

r x{t)dt. (8) 

^2 — JTi 

In general we will choose an initial time Ti large enough 
to eliminate transient behavior, and the interval T 2 — Ti 
large enough to reduce fluctuations. In order that the 
order parameters be small, we require that J2m=i ^nm ^ 
7V^. In particular, in our examples we have = Anm 
and so this condition becomes dn Under this 

assumption Rn = 0 for all n is an approximate solution 
of the system. In this section, we will study the stability 
of this incoherent state using a method similar to that 
used in Refs. [15, 16]. 

A. Dispersion relation and onset of instability 

With Rn = 0, (5) implies that Pn = 0 or Pn = pn = 
constant and each oscillator rotates with a constant an¬ 
gular frequency. Let us denote these solutions by 

9n{t)=Pnt + 9l , 

Pn{t) = Pn ■ 

We will assume that the initial phases are uniformly 
distributed in [0,27r). Letting {69n,6pn) denote small 
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Figure 
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5 
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5000 
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TABLE I: Parameters for the Erdos-Renyi (ER) networks studied here. Here N is the number of nodes, ||(i|| is the mean degree, 
and A is the largest eigenvalue of A. In all cases, the correlation coefficient, (17), is p = 1. 
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TABLE II: Parameters for the scale-free (SF) networks studied here. Here a and dmin are parameters in (16), A is the largest 
eigenvalue of A, and p is the correlation coefficient (17). In all cases the number of oscillators is = 10^ and the mean degree 
is ||ci|| = 100. 


perturbations to the incoherent state (0„(t),pri(t)), lin¬ 
earizing (2)-(3) gives 
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mean of a product of functions of independent variables 
is the product of their means, we can approximate (10) 
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upon neglecting {K/N)Y^^Anm cos 

m ^n) — 

0{'/d/N). These equations can be solved for the fastest 
growing mode using the new variables 


This is an eigenvalue equation; indeed, suppose that A is 
an eigenvalue of A and {&«} its corresponding eigenvec¬ 
tor, then (11) gives 
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As we show in App. A, upon setting Bn{t) = with 

the complex growth rate s = j + iuj, and assuming that 
7 > 0, then in the limit t > oo the eigenvector {t>„} and 
growth rate are determined by the eigenvalue problem 


K Aknbn 


( 10 ) 


Equivalently, 2N/K is an eigenvalue of the matrix 
Adiag{(s — ipn)~^}- For a given matrix A, distribution 
of initial momenta pn, and coupling constant K, (10) de¬ 
termines the growth rate 7 and oscillation frequency w 
of perturbations from the incoherent state. 

In the rest of the paper, we will consider—for 
simplicity—the case in which the initial momenta Pn are 
independent of the network properties, i.e., of A. That 
is, we can consider the set {{AkmPn) '■ n = 
to be a sample from a joint distribution that is, in fact, 
a product of two independent distributions, one for the 
network (A) and one for the initial conditions (p). In this 
case, we propose to look for solutions {bn} of (10) that 
are also statistically independent of the momenta. This 
hypothesis will be verified, a posteriori, below. Since the 


This verifies the hypothesis: if the momenta are uncor¬ 
related with A, they will also be uncorrelated with its 
eigenvectors {bn}, thus justifying our derivation of (11). 
The eigenvector that corresponds to the earliest onset of 
instability (i.e., the smallest K) is that corresponding to 
the eigenvalue of A with largest magnitude, which we 
will henceforth denote just by A (we assume A is non¬ 
negative, irreducible and aperiodic so that A is unique 
by the Perron-Frobenius theorem). In what follows we 
study the growth rate associated with this mode. If the 
momenta p have the distribution g{p) we can write, in 
the limit A —>■ 00, 

g{p)dp 

K\ y_oo (s - ipY 

Integrating by parts, 

r 9'{p)dp 
KX s-ip 

where g' = dg/dp. Now let s = 'y + ioj. Inserting this and 
separating real and imaginary parts and noting that A is 
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FIG. 1: Time average of the order parameter R, i.e., {R)t as a function of K [(a) and (c)] and XK [(b) and (d)] for a variety 
of network structures having N = 10"^. Panels (a) and (b) show the results for simulated Erdos-Renyi networks with varying 
q. Panels (c) and (d) show the results for simulated scale-free networks with varying a. From (fe) and (d), we see that plotting 
(R)t against XK for any network causes the transitions to line up, in agreement with (15). 


real since A is symmetric, we get: 

^ ^ r°° g'{p){uj-p)dp 

KX y_oo 72-b(w-p)2 
/■°° g'ip)dp 
7_oo l^ + {0J-p)^ 
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As an example, we consider the case in which g(p) is 
a Gaussian centered at p = 17 with standard deviation 
ctq. By symmetry, the second equation is satisfied when 
uj = n. The first equation in (13) then yields 
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We find the critical coupling strength Kc by letting 
7 —>■ O’*", obtaining 


Kc 


2alN 
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(15) 


The dependence of Kc on the largest eigenvalue of A is 
similar to that observed in various other dynamical sys¬ 
tems on networks such as the Kuramoto model [15], epi¬ 
demic spreading [17], and the propagation of avalanches 


[18]. The largest eigenvalue captures various effects of 
network structure including the degree distribution and 
degree-degree correlations [19]. 

To get the growth rate 7 for any given K > Kc, we can 
invert equation (14) numerically. We note that, given a 
value of K/Kc, the growth rate 7 is independent of the 
structure of network. 


B. Numerical experiments 

Here we present computations for several examples to 
validate and illustrate our results and verify some of the 
assumptions in the derivations of § 111 A for the linear 
stability of the incoherent state. As network examples, 
we use both Erdos-Renyi and scale-free networks with N 
nodes. We generate the Erdos-Renyi networks by estab¬ 
lishing an undirected link between nodes n and m (i.e., 
setting Anm = 1 ) with probability q, and not establishing 
a link {Anm = 0 ) with probability 1 — g to obtain a net¬ 
work with mean degree ||d|| = {N — l)q. The scale-free 
networks are generated, using the algorithm of Chung 
and Lu [20], to have a target degree distribution of the 
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FIG. 2: Time averaged global order parameter i? as a function of K and XK for scale-free networks with varying edge-degree 
correlation p. The three networks have distinct A [the inset of panel (b)], and panel (a) shows they have distinct Kc- Panel (b) 
plots {R}t against XK showing that the three curves nearly coincide as they did in Fig. 1. 
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FIG. 3: Time-averaged global order parameter for a fixed Erdos-Renyi network as a function of (a) K and (b) K/oq. The 
initial momentum distribution is a Gaussian with mean fl = 6 and varying standard deviation cro, as shown in the inset of (b). 


form 


P(d) 


d (i > ^min ; 

0 d djniri ■ 


(16) 


Given a value of a, we choose dmin to achieve a desired 
mean degree ||d||, which will be noted for each specific 
case. Tables I and II show the parameters used in the 
various experiments for the two types of networks. 

In the following set of experiments, we fix iV = 10^ for 
both network types. We start with the phases 0„(O) uni¬ 
formly distributed in [0, 27r). The initial momenta p„(0) 
are sampled from a Gaussian distribution with mean 
12 = 6 and a standard deviation ao = 0.35. We integrate 
Eqs. (2) and (5) using a second-order leap-frog algorithm 
with time step h = 0.01. In most of our experiments we 
report a time averaged value of the global order param¬ 
eter, i.e., {R)t as a function of the coupling strength, K. 
Integrations start at K = 0, and K is periodically in¬ 
cremented by AK (in the plots, AK is the separation 
between consecutive symbols). Integration at the new K 


value continues from the current state. The total inte¬ 
gration time for each value of K is typically T 2 = 1000 
time units; this includes a period during which transients 
decay, typically Ti = 500 time units. The time average 
of R is computed using (8). 

The first experiment studies the effect of varying the 
link probability q for the Erdos-Renyi networks and of 
varying the degree exponent a (16) for the scale-free net¬ 
works; results are shown in Fig. 1. Panel (a) shows the 
time averaged global order parameter (6) as a function of 
the coupling strength K for Erdos-Renyi networks with 
various values of q [indicated in the inset of Fig. 1(b)], and 
(c) shows the same quantity for scale-free networks with 
various values of a [indicated in the inset of Fig. 1(d)] 
with ||(i|| = 100. These networks have different largest 
eigenvalues, A (recall This. TII), and therefore, in agree¬ 
ment with (15), synchronization begins at different val¬ 
ues of K. However, for all of the networks the onset of 
the transition to synchrony begins at the same value of 
XK, see panels (5) and (d), confirming that Kc oc 1/A. 
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FIG. 4: Time averaged global order parameter as a function of (a) network size, N, and (b) edge probability, q for Erdos-Renyi 
networks with K = |ifc. The dashed lines have a slope of — | and an arbitrary intercept, corresponding to the theoretical 
estimate {R)t ~ {qN)~^^^. 


Remarkably, the entire set of curves collapse onto a sin¬ 
gle curve, indicating that even the partially synchronized 
states depend on the network structure only through A. 

As a second experiment, we use scale-free networks to 
study the effect of increasing degree-degree correlations. 
Correlations between the degrees dn, dm of the nodes 
connected by a randomly chosen edge, Anm ^ 0 , (also 
known as assortative mixing by degree [ 21 ]) can modify 
the eigenvalue A which, by (15), should affect the onset 
of synchrony. These correlations can be quantified by the 
coefficient [19] 

_ II dn dm lie 

||d„||e||dm||e ’ ^ ^ 

where || • ||e denotes an average over the edges. In our 
experiments we first construct a scale-free network with 
dmin = 33 and a = 2.5; then we rewire the edges fol¬ 
lowing the algorithm of [19] to increase p. The initial 
network has p = 0.78; subsequent rewiring creates an in¬ 
termediate network with p = 0.94, and a final network 
with p = 1.09. These three networks also have different 
principal eigenvalues, recall Tbl. II; however, all have the 
same degree distribution. The results, in Fig. 2 (a), show 
the onset of synchronization occurs at different values of 
AT, as in the first experiment. Therefore the simple scal¬ 
ing Kc oc l/|jd|| is not sufficient to describe the behavior 
of heterogeneous networks like scale-free networks. How¬ 
ever, when {R)t is plotted against ART, as shown in Fig. 2 
(b), the transition points again align as predicted by (15). 

In a third numerical experiment we study the effect of 
varying the distribution of initial momenta Pn- More 
specifically, using a single Erdos-Renyi network with 
N = 10^ and q = 0.01, we consider a Gaussian distri¬ 
bution of momenta g{p) with mean H = 6 and various 
standard deviations, aq. Figure 3(a) shows a plot of {R)t 
versus K for the different values of aq and panel (b) shows 
the collapsed version when the abscissa is Kld^. As ex¬ 
pected from (15), the critical values collapse to one point 


in the latter case, and—as before—the entire set of curves 
nearly coincide near the transition. 

As mentioned before, in the incoherent state we expect 
that Rn has fluctuations of order \/dn/N. Therefore R = 
Y^n^i^ri/WdW should scale as i? ~ |j^/^||/||c?||. For Erdos- 
Renyi networks the distribution is sharply peaked about 
||c?|| « Nq. Therefore we should have that R ~ {qN)~^^^. 
To verify that the observed finite value of R is consistent 
with these finite-size effects, we varied the network size 
N and edge probability q, holding the coupling constant 
fixed at AT = ^Kc- The results are shown in Fig. 4. The 
time average {R)t is shown as a function of N for q = 0.01 
in panel (a), and as function of q for N = 10"^ in panel (b). 
These results show that {R)t ~ and {R)t ~ 

respectively (the dashed lines have a slope of — | with 
arbitrary intercept), consistent with fluctuations due to 
the finite size of the network. 

Finally, we test our results for the linear growth rate 
7 from (14). Here we use an Erdos-Renyi network with 
N = 25,000 and ||d|| = 5000 (the reason for the larger 
N and \\d\\ is discussed below). We plot R{t) on a log 
scale as a function of t in Fig. 5. The solid lines represent 
data from direct numerical integration of (2)-(3), and the 
dashed lines have the slope 7 predicted from (14). 

Of course, exponential growth predicted by the lin¬ 
ear theory can occur only when R 1. In addition, 
(14) gives the growth rate of the fastest growing mode, 
but initial conditions may contain a mixture of different 
modes. Therefore, we expect the theoretical growth rate 
only over an intermediate time domain where the fastest 
growing mode dominates, but where R{t) has not yet sat¬ 
urated. To make his region as large as possible within our 
computational constraints we chose ||d|| = 5000. Since a 
quantitative comparison would require to arbitrarily se¬ 
lect an interval [AmiNjAmax] to compute a slope, we 
present here just the curves in Fig. 5 and do not attempt 
to fit portions of these curves. Despite these difficulties, 
the simulations of the full model show a growth rate that 
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FIG. 5: We study the growth of the order parameter i? on a 
log scale over time plotted on the a:-axis. The solid lines show 
simulated data for an Erdos-Renyi network and the dashed 
lines show the corresponding theoretical result from (14). The 
different curves are for a variety of K values. Theory and 
simulations agree well with increasing KjKc values. 


seems to be well approximated by (14). 


IV. SYNCHRONIZED STATE 

In §III, we studied the incoherent state; now we turn 
our attention to synchronized solutions, i.e., solutions for 
which the local order parameters R„ are nonzero even in 
the IV —>■ 00 limit. We are interested in the long-time av¬ 
erage of the order parameters (4), {Rn)t- Following [10] 
and based on numerical experiments (see below), we look 
for solutions such that the different local order parame¬ 
ters exhibit, on average, a phase synchrony—they rotate 
with a common frequency 14 and at a common phase (p: 

ItflO ~ / 72^ 5 

where is constant. This is a nontrivial assumption that 
we expect to be valid when all the nodes have neighbors 
that are representative of the network as a whole and 
the distribution of momenta is sufficiently narrow. To 
implement this assumption, we write 

= (r„ + , (18) 

where the real average local order parameters r„ and 
global phase ip are defined by 

(19) 

so that {Zn)t = 0. There are two implicit assumptions 
here: first, there is an 14 such that \{zn{t))t\ ^ r„, so 
that Zn represents the fluctuations, and second, there is 
a single phase p that makes all of the Vn real. The goal 
of this section is to obtain equations that will determine 
both the common frequency 14 and the local order pa¬ 
rameters Tn- 


As in [10], we define new variables and in a ro¬ 
tating frame, 


dn = 8n- (144 + p) 

Pn=Pn-^ ■ 


Inserting these and (18) into (2) and (5) gives 


On=Pn , (21) 

Pn = -Ktu sin(0„) -I- ATm( 2 ;„e“*®") . (22) 

The second term of (22) can be thought of as a perturba¬ 
tion to the Hamiltonian dynamics of each oscillator that 
preserves the total energy of all the oscillators (1). We 
treat this perturbation as if it were stochastic and assume 
that the probability of observing node n in a given region 
of the phase space (0„,p„) over a long time is given by a 
Boltzmann distribution [22]. More precisely, we assume 
that for any function / of the single oscillator variables, 
for large starting time Ti and large interval T 2 — T), the 
time average (8) limits to a phase space average: 

n2'K nOO 

{f)t {f)g = If {On,Pn) 9 (ffn,Pn',rn) dp„ d0„ 

Jo J — oo 

(23) 


Here g is the Boltzmann distribution for the single-rotor 
energy 


9(9,p;r) 


/ 31/2 

(27r)3/2Io(iV/3r) 


g-P{p^i'^-Kr cos(e)) 


(24) 


for an inverse temperature j3 that must be determined. 
The Bessel function, Iq, in the denominator normalizes 
the distribution: fQ^Jf°^ 9 (^,P',r) dpdd = 1. For this 
distribution, the mean square momentum (in this case, 
the variance of p) is 


r2n 




P^ 9{P, r) dB dp = /3' 


(25) 


> —OO 0 


and the mean potential energy is proportional to 



cos(0) g{0, p; r) d0 dp 


v{Kl3r) , 


(26) 


where we introduce the notation 


v{x) 


Ii(a^) 

Io(a:) 


(27) 


and Ii is the first order Bessel function. 

Using (19) and (20) in the definition (4) of the local 
order parameter, we can solve for r„, and then use (23) 
and (26) to obtain 


Tn 


1 

N 


N 

^ ^ -^n7n{^Os{9jri)}t 
m—1 


1 J 

AnmV{KPr) . 

m—1 

( 28 ) 








FIG. 6: Time averaged global order parameter, {R)t as a function of the coupling strength K for (a) Erdos-Renyi {q — 0.01) 
and (b) scale-free {a = 2.5, dmin = 33) networks with N = 10^ and ||d|j = 100. The circles (blue) are the result of numerical 
integration for an initial Gaussian distribution go{p) with mean Pq = 6 and variance ctq = 0.12 (for Erdos-Renyi) and ctq = 1 
(for scale-free) and random phases uniformly distributed in [0, 27r). The solid curves (black) are the numerical solutions to (33). 
The dashed curves (red) show the second-order approximation (38), valid near K = Kc (see §IVA). 


Equation (28) depends on the inverse temperature (3 
introduced in (24), which can be determined by conser¬ 
vation of energy. Suppose that initially the rotors have 
a distribution of momenta with mean Pq = |jpra(0)|| and 
variance CTo = ll(Pn(0) ~ -Po)^||, and that they have a 
distribution of phases 0„(O) with potential energy Uq = 
-Kj{2N) J2n,m cos(0„(O) - 6»„(0)). The initial en¬ 

ergy is then 

Eq = — {Pq + o-q) + Uq . (29) 

Since the total momentum is conserved by (3), the mean 
momentum at any time remains equal to Pq. Under the 
Boltzmann assumption (23), the mean {p)g = 0, which, 
by (20), implies that 

Po = {Q+p)g = n . (30) 

In the new coordinates (20), the total energy (1) at 
any time is 

N N N 

E = ^ ^ (U -b p„)^-y ^ r„ cos (0„)- Y ^ Re(z„e"*® 

n—1 n—1 n—1 

Since the energy is constant, we can take a time average 
and use (25) to obtain 

N K . - K . 

E = y (u^-bcr^)-— r„(cos (6»„))t-y Y 

n—1 n—1 

(31) 

We now neglect the terms proportional to the fluctua¬ 
tions Zn (see below for a discussion). Since energy and 
momentum are conserved, Eq = E and Pq = fl, we can 
apply the Boltzmann assumption (23) and combine (26), 
(29), and (31) to compute the variance: 

K . 2 

= ctq + — ^ r^v{Kl3r) + —Uq . (32) 


Substituting for (3 using (25) in (28) and (32) gives a 
closed system of iV -|- 1 equations for the local order pa¬ 
rameters and the variance: 



This system generalizes analogous self-consistent re¬ 
sults for the all-to-all coupled case (e.g., see Eq. (16) in 
[4]) Note that this system always has the trivial, incoher¬ 
ent solution r„ = 0, n = 1,..., N' and = (Jq + 2Uq/N. 
By the analysis of §IIIA, this solution is stable when 
K < Kc- We note that when the initial conditions are in 
the incoherent state, i.e., when the phases are uniformly 
distributed in [0, 27r), the potential energy term 2Uq/N 
is negligible in the limit N ^ oo. 

) ^ To find a nontrivial, synchronized solution with > 0, 
we solve the system (33) numerically for the iV -|-1 vari¬ 
ables {rn} and cr^. A simple method is fixed point iter¬ 
ation: given a guess {r„ > 0} and cr^, new values can be 
computed from (33) using the guesses on the right-hand 
sides. Numerically, this iteration converges to values that 
t appear to be independent of the initial guess, suggesting 
that there is a unique solution to these equations, and 
that there is a nontrivial solution, r„ > 0, when K > Kc- 
Once the local order parameters are known, the global 
order parameter is computed from (6). 

In Eig. 6 we show a comparison of the predictions of 
(33) (black solid lines) with direct numerical integration 
of the ODEs (2)-(3) (blue circles) for both an Erdos- 
Renyi network, panel (a), and a scale-free network, panel 
(b). The theory agrees with the simulations, except that 
when K < Kc the observed order parameter R is not 
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FIG. 7: Time averaged oscillator frequencies, panel (a), and fluctuation amplitudes, panel (b), for oscillators in an Erdos-Renyi 
network (q = 0.01, ||d|j = 100, N = 10"^), with K = 600 = 3Kc. The initial momenta had a Gaussian distribution with mean 
Pq = 6 , and variance (Tq = 1. Time averages were taken over an interval [2000,10000]. Panel (b) shows the time average of the 
fluctuating terms in (31) relative to the synchronized terms. 


zero, as predicted by the theory, because of finite-size 
effects. 

Before moving on to the next section, we present a 
discussion of the three main assumptions made to derive 
(33): 

Asmp. 1. The local order parameters have a eommon 
rotation frequency and phase, t/jn = ^t + ip [introdueed in 
(18)]. 

Asmp. 2. The final state is ergodic and has a Boltzmann 
distribution (23) [used in (28) and (31)]. 

Asmp. 3. Fluctuations can be neglected: 
|j(Re( 2 ;„e"*®"))t|| < || (?-„ cos (0„))t|| [for (31)]. 

The first assumption, that the phases of the local order 
parameters are all the same, is reasonable when the net¬ 
work is constructed in such a way that the neighbors of 



'tpnit) - l(l{t) 


FIG. 8: Histogram of'i/)„(t)— i/)i (t), for the network of Fig. 7a, 
but with N = 10^ oscillators. Here the are sampled each 
unit of time over a total of 2000 time units, so the total num¬ 
ber of events is 2(10)®. 


different nodes have the same statistical properties and 
the initial momentum distribution is sufficiently narrow 
(since On = Pn). Indeed, this assumption has also been 
used successfully in studies of Kuramoto oscillators on 
complex networks [15, 23]. It is expected to break down 
for networks where the oscillator properties are correlated 
with the network structure, such as lattices with spatially 
dependent frequencies [24] and communities with differ¬ 
ent oscillator properties [23, 25-27], or when the distri¬ 
bution of momenta is bimodal [10]. 

Both the Erdos-Renyi and scale-free networks sat¬ 
isfy the statistical equivalence property. The validity of 
Asmp. 1 can verified numerically. For each rotor n, we 
can estimate its effective angular velocity by a time 
average, i.e., we compute 

= rp ^ rp bPn{T2) - 

for large Ti and T 2 — T 1 . As usual, Ti is chosen to elimi¬ 
nate initial transients and T 2 to decrease the noise. Typ¬ 
ical values are Ti = 2000 and T 2 = 10^. An illustration 
for an Erdos-Renyi network is shown in Fig. 7a. The 
figure shows that the deviations of Tin from the aver¬ 
age TI are of order 0.5% (and they become smaller as T 2 
is made larger). Even if the rotors have the same fre¬ 
quency, they could have different phases. To verify this 
is not the case, Fig. 8 shows a histogram of tpn {t) — '0i (t ), 
where the histogram samples all n yf 1 and at the integer 
times t = 1,2,..., 2000. The plot shows that the phases 
remain very close to each other. The tails of the distribu¬ 
tion correspond to phase slips [i.e., ipnit) — ipi(t) rapidly 
changing by 27r]; these slips become less frequent as the 
mean degree of the network is increased (not shown). 

Assumption 2 is commonly used in the analysis of the 
HMF model [1]. The idea is that the single rotor, de¬ 
scribed by (21)-(22), is a Hamiltonian system exchanging 
energy with the rest of the network, which for large N 
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FIG. 9: Single oscillator phase space distribution for an Erdos-Renyi network (A^ = 100, |jd|| = 10) with K = 1000 = 50Kc 
with an initial Gaussian distribution of momenta with mean Po = 6 and variance (Jq = 1 and uniform initial phases. Panel 
(a) shows the theoretical Boltzmann distribution (24) for node n = 1, with (3 = 0.965 and r\ = 0.036. Panel (b) shows the 
numerical distribution of averaged over the time interval [2000,10000]. The histograms in both panels use bins 

of AS = 0.0628 and Ap = 0.0718. 




can be taken to be a thermal bath. This implies that, 
in equilibrium, the statistical behavior can be described 
by the Boltzmann distribution (24). The validity of this 
assumption is demonstrated in Fig. 9, which compares 
the theoretical distribution g(0, p;r„) with calculated 
from (33) in panel (a), with a histogram of the empirical 
long-term distribution of the variables (0„,p„), in panel 
(b) for an Erdos-Renyi network. The figure shows the 
results for the arbitrarily chosen node n = 1—results are 
similar for other choices. For the phase space of the cho¬ 
sen rotor, the theoretical and experimental distributions 
are visually close. 

Assumption 3 would follow if the fluctuations were 
uncorrelated with 0„, because then = 0, since 

the fluctuations satisfy, by definition, {zn{t))t = 0. If 
the number of connections per node is large these cor¬ 
relations should be weak, since Zn is determined by the 
behavior of Om for all the neighbors m of node n, each of 
which, in turn, depends on the phases of all of their many 
neighbors. These heuristic arguments can be validated 
numerically by computing explicitly (Re( 2 :„e“*®"))t. In 
Fig. 7b we plot the ratio Re(z„e“*®'‘)t/|| (r„ cos(0„))t|| for 
each of the N = 10^ nodes in an Erdos-Renyi network. 
For most nodes this ratio is small, less than 0.1, though 
for 20 nodes it is larger than 0.1 and the maximum ratio 
is 0.16. The validity of Asmp. 3 depends on the network 
average of the numerator being relatively small, and for 
this case we found 


||Re(z„e-*®")t|| = l.l(10)-5 < ||(r„ cos(0„))t|| = 2.3(10)-3, 


confirming the assumption. 


A. Perturbative approximation 

Equations (33) allow us to calculate the order param¬ 
eter R given a network adjacency matrix A, a coupling 
strength K, and the total energy. Though the numeri¬ 
cal solutions for and a agree well with the simulations, 
they do not offer additional insights into how the network 
structure affects the general properties of the transition 
to synchrony. In this section, we present a perturbative 
analysis of (33) near the bifurcation at K = that 
allows us to determine what properties of the network 
affect the value of the order parameter close to the bifur¬ 
cation. 

In order to do this, we will solve this system perturba- 
tively, assuming that AAT = AT — ATc <C 1—the coupling 
constant is just beyond the critical value. We solve the 
system (33), i.e.. 



using ^ for the variance and setting the initial po¬ 
tential energy Uq to zero, e.g., for initial conditions in 
the incoherent state in the limit A —>■ oo. Introducing 
a formal small parameter e, the perturbative expansion 
takes the form 

K = Kc + eAK, 

Tn = + C’(e^), (34) 

/r = -I- 

where we have anticipated already that r„ ~ (AA")^/^ 
and have included only terms up to the order necessary 
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to determine Vn^ in the analysis that follows. Inserting 
these in (33) and expanding in powers of e we obtain at 
zeroth order, 






2 

0> 


as expected. 


The next terms, of order imply 



Kc 

2alN^ 


r(i) = 


0 , 


which gives 


K, 


2alN 

A 


= Cu. 


(35) 


Here u and A are the principal eigenvector and eigenvalue 
of A, and C is a normalization constant to be determined 
(as we will see, the product Cu does not depend on the 
normalization of u) . This result is in agreement with the 
linear stability calculation of §III [cf. (15)]. The terms 
of order lead to oc u (although this will not be 
used), and to 


^ ^|j„2|| ^ (3g) 

ZO-Q 

Here denotes the vector with components and the 
|j...|| the network average (7), as usual. Finally, the terms 
of order yield 



K, 

2alN 


A 


r(3) = 


C 

IGaljN 


A 


8(Tl{AKal - 




(37) 


In order to eliminate the unknown vector we mul¬ 
tiply (37) on the left by Since = A, then 

u^A = Au^ and using (35), the left hand side vanishes, 
giving the solvability condition 

0 = 8al{AKal - 


Using (36), this determines C: 

f sap^iiu^ii 
'vllu4|H-4|lu2p KiJ 

Finally, using the definition (6), 




and (35), we find (dropping the formal parameter e), the 
main result of this section. 



A / 2||u^||||uf 

- IMII 'vllu4||+4||u2p 


1/2 


(38) 


This expression provides some insight into the effect 
of heterogeneity on synchronization through the factor 
G. For example, for an uncorrelated network for which 
Un OC dn [19], G = y/2/5 for a regular, homogeneous 
graph with dn = d, while G —>■ 0 for when the degree dis¬ 
tribution is heterogeneous so that Hd'^ll — >■ oo in the limit 
TV — >■ oo. Thus, in this case we find that heterogeneity 
tends to make the transition to synchrony less sharp. 

To illustrate this, compare this theoretical result to the 
numerical results for the time averaged order parameter 
{R)t in Fig. 6 for an Erdos-Renyi (homogeneous) and 
a scale-free (heterogeneous) network with the same size 
and mean degree. The dashed lines show the approx¬ 
imation (38). We find G = 0.633 for the Erdos-Renyi 
network, larger than G = 0.354 for the scale-free net¬ 
work, as we would expect. The second-order approxima¬ 
tion agrees with the numerical results for the scale-free 
network whenever K > K^', however, R for the Erdos- 
Renyi network is not well approximated for higher values 
of AK. Of course, AK is assumed to be small in the 
derivation above, so there is no reason for agreement for 
large AK. 


B. Large K limit 


Figures 3 and 7 suggest that R tends to an asymptotic 
value, R—>'R<I, asiV—J-oo. In this section we 
will study this limit and explore how R depends on the 
network. To begin our analysis, we divide (33) by K and 
let r] = A/K to obtain 



Under the hypothesis that lim^_,.oo i] = fj is finite, in the 
limit K ^ oo the system above reduces to 


r 


n 


1 X V 

(/V X! 




n—1 



(39) 


for the asymptotic values f„, and t). A solution can be 
found numerically as the fixed point of the relaxed iter¬ 
ation scheme 


771^1 ^ ^ 


R« G 















12 


1 

0.9 
0.8 
0.7 
« 0.6 
0.5 
"^ 0.4 
0.3 
0.2 
0.1 
0 

2.5 3 3.5 4 

a 

FIG. 10: Asymptotic order parameter, R, for scale-free net¬ 
works with varying degree exponent a. The red triangles show 
the iterative solution of (39), and the blue circles show {R)t 
from simulations using a large coupling constant, K = 10^. 
The dashed line shows the approximation from (42). 


where 0 < ,5 < 1 is a relaxation factor included to obtain 
convergence to a fixed point. 

Figure 10 shows the global order parameter R obtained 
from both the numerical solution of (39) (red triangles) 
and the numerical solution of the full system, (2)-(3) 
(blue circles) for scale-free networks as a function of the 
exponent a of the degree distribution (16). For the sim¬ 
ulations, K = 10^ and the observed order parameters 
are about 1% smaller than the theory, which assumes 
N ^ oo This difference is not visible on the scale of the 
figure. Remarkably, the asymptotic value of the order 
parameter is nearly independent of the exponent a. 

An additional approximation can be made, following 
[28], if we assume that the local order parameters are 
proportional to the nodal degrees, = Bdn, where B is 
a constant to be determined. This approximation works 
well for homogeneous networks without correlations [19]; 
for scale-free networks, it works well for power law expo¬ 
nents a > 3 [15]. If we replace by Bdn in (39) and sum 
the first equation over n we obtain, using d^ = ^nm, 
that 

= 4 ^ (40) 

rr, = 1 \ 7 / 



Comparing the two equations we see that rj = B^N\\d\\. 
From the definition of the global order parameter we find 
R = J2n^n/\\d\\ = BJ2ndn/\\d\\ = NB. Using these two 
in (40) we find a nonlinear equation for the single variable 
R 

i?||d|| = — y . (42) 

UiMiiy 


This equation can be solved numerically using standard 
root-finding tools, and produces the dashed line shown 
in Fig. 10. 

If the degrees of individual nodes are not known, but 
the degree distribution P{k) is known, one can approxi¬ 
mate the sum in (42) by an integral to obtain an implicit 
equation for R (using the dummy variable k instead of 
d) 

Finally we note that, although this was not pursued here, 
a similar mean mean field approach (i.e., r„ = Bdn) could 
be used to further study the system (33). 

V. CONCLUSION 

In this paper we studied the HMF model where the 
interactions between rotors are described by a weighted 
adjacency matrix. We found that, as in other dynamical 
systems on networks (e.g., [15], [17], [18]), the transition 
to synchrony occurs at a value of the coupling constant 
inversely proportional to the largest eigenvalue. A, of the 
adjacency matrix. Thus the primary effect of network 
structure on this aspect of the dynamics is A. 

We obtained a set of equations that determine the 
set of local order parameters in the synchronized state. 
These equations relied on three assumptions that were 
verified a posteriori for the Erd os-Renyi and scale free 
networks studied in Sec. IV. Of these assumptions, the 
most important is that the network is constructed in such 
a way that the neighbors of all nodes share the same 
statistical properties. This assumption is not satisfied, 
for example, by networks with strong community struc¬ 
ture. While this seems restrictive, the class of networks 
to which our results apply include networks with hetero¬ 
geneous degree distributions (e.g., scale-free networks) 
and networks with degree-degree correlations. It is also 
expected that some of our results could be extended to 
networks with community structure. 

Our main result is a method to quantitatively explore 
the effect of network heterogeneity on the transition to 
synchrony, resulting in (38). In addition to determining 
Kc, the critical coupling constant, network heterogeneity 
also affects the sharpness of the transition, with more het¬ 
erogeneous networks having a less pronounced transition. 
However, even though such heterogeneity (as represented 
the degree distribution) has a strong effect both on the 
location and sharpness of the onset of synchrony, it seems 
to have little effect on the degree of synchronization for 
large coupling. 

In conclusion, our results show that many of the phe¬ 
nomena that have been observed for the all-to-all coupled 
HMF model persist for more complex networks, and that 
the onset of synchrony is determined by spectral proper¬ 
ties of the coupling matrix. 
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To find the dispersion relation, we assume exponential 
growth of the perturbations, i.e., Bk{t) = where 

Re(s) > 0. Using this, we get that 


Appendix A: Derivation of the dispersion relation 


In this Appendix we derive the dispersion relation (10) 
by studying the evolution of the perturbations (SOn, Spn) 
to the incoherent initial state {0n{t),pn{t)) = {pnt + 
), where are uniformly distributed in [0,27r), and 
the initial momenta are arbitrary. Integrating the per¬ 
turbed ODEs (9) formally with respect to time gives 


59n{t) 

Spn{t) 


= / Spnit')dt' + S9nito) , 


'to 


K 


N 


^ Anm / COs[6»,„(f') - 9n{t')]S9m{t')dt' 

m=l ”'*0 

Spnito) , 


Defining C„(f) = 59n{t) — S9n{to) and integrating the 
second equation from Iq to t, we find 


n=l “'*0 Jto 

n—1 

(A2) 

Since we are assuming Re(s) > 0, the left hand side of 
(A2) grows exponentially with t. However, the term / 
grows at most quadratically. 


|/| < to f^ X! ^nmS9m{to) + {t- to)6pn{to), 


N 


N 


CM 


rt rt' N 

c / ^ -c^nm [ 

'to Jto 




(Al) 




and therefore as < —>■ oo the first term on the right-hand 
side of (A2) must balance the left hand side. Replacing 
9n = Pnt + 9n(0), we obtain 


where / is 


K ^ f* 

I Anm / / COs[9mit'') - 9n{t”)] S9m{to)dt”dt' 

m=l "'‘o Jto 

+ {t - to)6pr,{to) . 

Defining Rfe(t) = multiplying 

equation (Al) by and summing over n yields 


K 


N 


“ 27 V ^knbn 


t rt' 


(A3) 


to Jto 


K 

m 


N 




2ie^{0) 


t rt' 


n—1 


to Jto 


(A4) 


pt rt' 

Bk{t) = ^^Aun / 

„=1 Jto Jto 

N 


Since the angles 6>n(0) are uniformly distributed in [0,27r), 
the second term can be neglected for large N. Integrating 
the first term and taking the limit t —t oo with Re(s) > 0 
we finally obtain (10). 
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